Complete Reconstruction of the Wavefunction of a Reacting Molecule 
by Four- Wave Mixing Spectroscopy 
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Probing the real time dynamics of a reacting molecule remains one of the central challenges in 
chemistry. In this letter we show how the time-dependent wavefunction of an excited-state reacting 
molecule can be completely reconstructed from resonant coherent anti-Stokes Raman spectroscopy. 
The method assumes knowledge of the ground-state potential but not of any excited-state potential, 
although we show that the latter can be computed once the time-dependent excited-state wavefunc- 
tion is known. The formulation applies to polyatomics as well as diatomics and to bound as well 
as dissociative excited potentials. We demonstrate the method on the Li2 molecule with its bound 
first excited-state, and on a model Li2-like system with a dissociative excited state potential. 

PACS numbers: 03.65.Wj, 31.50.Df, 82.53.-k, 78.47.nj 



For several decades now, femtosecond pump-probe 
spectroscopies have been employed to study transition 
states of molecules reacting on excited potential surfaces 
[ll-Q. Although these studies have shed a tremendous 
amount of light on excited-state dynamics, none of the 
methods in use provides complete information on the 
excited-state wavefunction. The need for an experimen- 
tal method that will provide this information is com- 
pounded by the fact that theoretical ab initio calculations 
for excited states are difficult and of limited accuracy. 

There have been several theoretical proposals for com- 
plete reconstruction of an excited-state molecular wave- 
function from spectroscopic signals These studies, 
however, generally assume that one or more excited-state 
potentials (or the corresponding vibrational eigenstates) 
is known. A notable exception is a recently developed 
iterative method for excited-state potential reconstruc- 
tion from electronic transition dipole matrix elements [8] 
but this method does not appear to be applicable to 
dissociative potentials. Experimental work has focused 
on wavepacket interferometry of vibrational wavepackets 
[^, [l3| as well as electronic Rydberg wavepackets jll|, ■ 

The approach we present here assumes knowledge of 
the ground-state potential but not of any excited poten- 
tial. In principle, the approach is completely general for 
polyatomics. Our strategy is to express the reacting- 
molecule wavefunction, \9{t))^ as a superposition of the 
vibrational eigenstates of the ground-state Hamiltonian, 



(1) 



Since the vibrational eigenstates {iV'g)} are assumed 
known, the challenge is to find the time-dependent su- 
perposition coefficients Cg{t). 

Consider a two-state molecular system within the 
Born-Oppenheimer approximation. The nuclear Hamil- 
tonians Hg and He correspond, respectively, to the 
(known) ground and (unknown) excited potentials, which 



can be of any dimension. For simplicity, we consider 
a 5-pulse excitation as well as a coordinate-independent 
electronic transition dipole, fi (Condon approximation). 
Applying first-order time-dependent perturbation theory, 
the wavepacket that we want to reconstruct is 



13| 



\^^>{t))^-^e-'"''{~^lel}\^Jo), 



(2) 



where the initial state, IV'o), is the vibrational ground- 
state of Hg with the eigenfrequency wo, £i is the ampli- 
tude of the pulse and t is the propagation time on the 
excited state measured from the time of pulse excitation. 
(Here and henceforth we take H = 1.) 

Substituting Eq. ^ into the definition of Cg{t), we 
find that the superposition coefficients are given by 



Cg{t) = i^ei(V'g|e ^-^^'IV-o) = il^£iCg{t) 



(3) 



Hence, the central quantities required for reconstructing 
5'(i)) are the cross-correlation functions Cg{t). It has 
long been recognized that these correlation functions ap- 
pear (up to /x) in the time-dependent formulation of reso- 
nance Raman scattering (RRS) [14i]; however, the exper- 
imental RRS signal involves the absolute-value-squared 
of the half- Fourier transform of the correlation function, 
hence the latter cannot be recovered from that signal. 

Fully resonant coherent anti-Stokes Raman scattering 
(CARS) has been shown to be a powerful probe of ground 
and excited electronic states properties [1^, [3| • In this 
letter we show that the correlation functions {cg{t)} 
may be completely recovered from femtosecond reso- 
nant CARS spectroscopy, allowing complete reconstruc- 
tion of the excited-state wavepacket. The formula for 
the CARS signal produced by a three-pulse pump-dump- 
pump sequence is P^^\t) = {'il)'^'^\T)\jl\ili^^\T)) + c.c. 
[ITj . where ^/^'^^^(r) is the third-order wavefunction and 
^(o)(r) = e~^^o'^ipo- Within the (5-pulse and Condon ap- 
proximations, P*-^-* takes the form 



-iH^T2 



'li^o), (4) 
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where r^j = r.^ — Tj is the (positive) time-delay between theorem we obtain a sine-type of spectrum with peaks at 



the centers of the ith and jth pulses and r43 = t — 
with r being the time of signal measurement. We have 
denoted Hg = Hg - ujq, e = z3^4^^^^^^gj(^o(^2i-HT43) ^j^h 
£1,2,3 as the first, second and third pulse amplitudes, re- 
spectively, and T = [t2i, T32, T43]. In writing P'^^^t) as a 
complex quantity we have assumed the signal is measured 
in a heterodyne fashion. 

As illustrated in Fig. [1] Eq. ^ has the following phys- 
ical interpretation: A first laser pulse, the pump pulse, 
transfers amplitude to the excited potential surface cre- 
ating a wavepacket whose time-dependence we are in- 
terested in reconstructing. After evolving on the excited 
state for some time, a second laser pulse, the dump pulse, 
transfers part of this amplitude back to the ground state 
where it evolves for a second interval of time. Finally, a 
third laser pulse excites part of the second-order ampli- 
tude to the excited state, generating the third-order po- 
larization that produces the CARS signal, measured at 
later times. The desired wavefunction |^'(i)) (Eq. ([2])) 
may already be recognized in Eq. (jlj. 




Molecular coordinate 

FIG. 1: Color online. The pump-dump-pump CARS scheme. 
^{t) is the desired wavefunction. 

The reconstruction of |^'(i)) from P'^^^ proceeds in five 
steps: 

1. Insert a complete set of ground vibrational states. 
Introducing J^g = 1 into Eq. (jH) we obtain 

N 

P(3)(T) =e^e-'=«^-Pf)(r43,r2i), (5) 

9=0 



(3) 

where Pg '(t43,T2i) = (V'o|( 



IV's)(V'ff|e' 



and ujg = (jjg ^ ujo- N is determined by the number 
of ground vibrational states required to expand |^'(t)). 
Note that the desired correlation functions Cg{t) may al- 

(3) 

ready be recognized in Pg . 

2. Fourier-transform P*^^^ with respect to T32. The 

transformation resolves P^'^^ into individual ground-state 
(3) 

components, Pg . Since T32 is defined to be positive, we 
multiply Eq. ([5]), prior to the transformation, by the 
rectangular function that takes the value 1 for the T32 
domain and elsewhere. Using the Fourier convolution 



the frequencies uj 



N 



P(''(T43,w,T2i) = ^5(a.,.g)Pf )(r43,T2i), (6) 

3=0 

where S{uj, g) = 2Tee''^'^-^o^'^^^^+^hmc[{uj - uJg)T], 
2T — T32 — T32, and f32 (T32) is the minimal (maximal) 
value of r32. Fixing (t43,T2i), Eq. ([5]) can be written as 
a matrix equation: 



p(3) _ o p(3) 



(7) 



3. Invert the matrix equation ([^. The equation 
Pg^' = Sgj^Pif' isolates the two-dimensional functions 

Pg (''43 7 ''"21) • In inverting S we choose the number of 
frequency elements (lu) equal to the number of the g el- 
ements so that the matrix is square. For numerical ac- 
curacy, the inversion is implemented separately around 
each of the peaks at ujg. 

(3) 

4. Take the square-root of Pg . Assuming the func- 
tions {"ipgix)} are real, we can rewrite Pg as 



PfHT43,r2l) 



(^,|e-^=--|V^o)(^3|e-^''=^^1V^o). (8) 

Taking the square-root of the diagonal of Pg {t43,T2i) 
(i.e. r43 — T21 — t), we recover the Cg{t) up to a sign: 



P^'\t) ^ ag{^g\e 



IV-o), (9) 



where Og = ±1 and the sign of ijjg{x) is as yet undeter- 
mined. By demanding continuity of the cross-correlation 
functions (and their derivatives), the coefficients Og can 
be regarded as time- independent. Substituting Eq. (|9]) 
instead of Cg (t) into Eq. ([3|) and using the resulting Cg (t) 
in Eq. ^ yields 
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\-^{t)) ^ifieiY,\^9)(^9\e~ 



1^0). 



(10) 



The different sign combinations of ipg{x) generate 2^+-'^ 
possible superpositions. (In fact, only 2^ are physically 
meaningful since we are free to set the sign of one of the 
f/-components.) Only one out of the 2^ |^(i)) coincides 
with |^(i)): the |^(i)) for which the sign combination 
satisfies Y.g = 1- 

5. Discriminating l'I'(t)) from the set {|5'(i))}. The 
set of wavefunctions {|^(t))} are all consistent with the 
CARS signal at a specific value of T43 = r2i [3]. How- 
ever, only one |4'(f)) is consistent with the signal deriva- 
tives. To see this, consider the rath derivative of the ex- 
perimental signal, Eq. with respect to T21: 

9"P(3)(r) 



:et(**(T43)|e-'^«^-i?^'|*(r2i)) 



9,9' 



-itAlgT32 



Cg{T43)Cg' {T2i)H"gg, , 



(11) 
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where 



T41 = r - Ti, 
= iH,-LOo)", and if^^^^, = (V'gliJrlV^g')- Substitut- 
ing instead of into Eq. ([TTl) gives 



= et 



(12) 



Accordingly, the I'^it)} for which ^"^l'y> = 

for all n, is the wavefunction that coincides with \'^{t)) 

of Eq. ([2]), and hence, is the reconstruction solution. 

In practice, we proceed as follows. We invert the time- 
dependent Schrodinger equation to calculate a set of po- 
tentials from each 



V{x) 



1 



d 



1 

2m dx^ 



(13) 



where m is the system's reduced mass. One can show 
that the potentials calculated by the |4'(i)) that do not 
coincide with \^{t)), are time-dependent fl9|. Only 
the potential calculated with |^(i)) = |4'(i)) is time- 
independent and hence corresponds to the excited-state 
Hamiltonian of the measured system. Thus, in order 
to find the correct wavefunction we use the set of cal- 
culated potentials, as if they were time- independent, to 
propagate the corresponding back to time zero. 

Of all the potentials, only the truly time- independent one 
will propagate the corresponding \^{t)) correctly back to 
I'i/'o), and therefore this |^(i)) is the correct wavefunc- 
tion. Note that the above procedure requires knowing 
the signal as a function only of T32 and T21 — T43. 

TABLE I: The parameters, in atomic units, for the X, A and A 
potentials used in simulating the CARS signals. 





X 


A 


A 


D 


0.0378492 


0.0426108 


9.11267x10"^ 


h 


0.4730844 


0.3175063 


1.5875317 




5.0493478 


5.8713786 


7.3699313 


T 





0.0640074 


0.0640074 



To test the above reconstruction methodology, 
we simulated the CARS signal by calculating 
('0(°)(T)|/t|?/''^'^)(r)) as a function of the three time- 
delays, for two one-dimensional systems. The first 
is the Jj\2 molecule, with its ground {X) and first- 
excited {A) electronic states as Morse-type potentials, 
V{x) = D{1 - e-''(^-^o))2 + T. The second system, 
henceforth denoted d-Li2, has the Li2 ground state 
(X) but a dissociative excited potential of the form 
V{x) = De~^^^~^°^ +T (denoted A). Table U gives the 
potential parameters in atomic units used for the simu- 
lations. The parameters for the Morse- type potentials 
are based on data published in 

The wavepacket propagations employed in simulating 
p(3) were performed using the split-operator method 21 1 



on a spatial grid of 256 points in the range of 2-12a.u. 
with time spacing of At=0.1fs. A constant transition- 
dipole of 2a. u. was used, and the pulse amplitudes £1,2,3 
were taken to be 10~^a.u. The ranges of time-delay for 
the Li2 (d-Li2) system were T2i,43 = — 200fs (0 — 80fs) 
with spacing of 0.2fs. For both systems, we took T32 = 
3 — 6000fs with Ifs spacing. 

For Li2, we inverted Eq. ([7]) for each of the first 25 
peaks of P^{ijj) using the matrices S with 25 frequency 
grid points centered around the peaks at u)g. This pro- 

f3) 

duced 25 two-dimensional functions Pg , g = 0, . . . , 24. 

For d-Li2 the procedure was performed for the first 

(3) 

40 peaks, producing 40 two-dimensional functions Pg , 
g = 0,...,39. 
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FIG. 2: Color online. Snapshots of the real part of the recon- 
structed (circles, red) vs. the exact (dots, blue) wavefunction, at 
various times on the excited (A) potential (solid line) of Li2 . 

In Figs. [Hand 13] we present snapshots of the real part 
of the reconstructed first-order wavefunction for the Li2 
and the d-Li2 molecules, respectively. For Li2 (d-Li2) we 
superpose the first 25 (40) eigenfunctions ipgi^) using the 
cross-correlation functions obtained by the CARS analy- 
sis and maintaining ^ = I- The reconstructed 
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FIG. 3: Color online. Snapshots of the real part of the recon- 
structed (circles, red) vs. the exact (dots, blue) wavefunction, at 
various times on the excited (A) potential (solid line) of d-Li2. 

wavefunctions are seen to be in excellent agreement with 
the exact ones, obtained by direct calculation of the first- 
order population, for all propagation times. For the Li2 
system, a high quality reconstruction is already obtained 
by superposing just 20 basis functions. 
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FIG. 4: Color online. The reconstructed (circles, red) vs. the 
exact (dots, blue) A potential of Li2. 
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FIG. 5: Color online. The reconstructed (circles, red) vs. the 
exact (dots, blue) A potential of d-Li2. 

Having determined the wavefunctions we calculate the 
corresponding excited potential surfaces from Eq. (fT^ 
using eight-point (three-point) central finite-differencing 
for the time (spatial) derivatives. The time-step used 
was 0.2fs but very good results were also obtained us- 
ing 0.5fs. Figures |4] and [5] compare the reconstructed vs. 
the exact potentials. The wavefunction (absolute value) 
used in calculating the potential is shown by a black solid 
line. Note from Figs. |4] and [5] that combining the re- 
constructed potential from two points in time (e.g. 5 
and 70fs for Li2 and 5 and 79fs for d-Li2) is sufficient to 
reconstruct the potential over the full range of interest 
(2-5A). Once the potential is known one can calculate 
the excited-state wavefunction as a function of time for 
any excitation pulse sequence without the need for any 
additional laboratory experiments. 

To conclude, we have presented a methodology for 
the complete reconstruction of the excited-state wave- 
function of a reacting molecule by analyzing a multi- 
dimensional resonant CARS signal. The methodol- 
ogy applies to polyatomics as well as diatomics. We 
have assumed that only the ground-state potential is 
known. The approach is very compelling since the de- 
sired excited-state wavefunction is explicitly contained 
in the formula for the CARS signal. Highly accurate 
reconstruction is obtained even far from the Franck- 
Condon region. In fact, in practice the method may 
be more accurate far from the Franck-Condon region. 



since the frequency shift between the pump and dump 
pulses will be more effective in discriminating unwanted 
processes that may contribute to the measured signal at 
k = ki — k2 -I- ka. We simplified matters by considering 
(5-function pulse excitations, a coordinate-independent 
transition dipole moment and only one excited-state po- 
tential. In future work we will test the removal of all 
these assumptions. 

We have shown that once the time-dependent wave- 
function is found, the excited potential can be recon- 
structed with quite high accuracy. It will be of great 
interest to test the method on polyatomics, where ob- 
taining multidimensional potential surfaces from spectro- 
scopic data has been one of the longstanding challenges 
of molecular spectroscopy. An important application of 
excited-state potential reconstruction will be the ab ini- 
tio simulations of laser control of chemical bond breaking. 
Experimental laser control has been greatly hindered by 
the lack of detailed theoretical guidance, which in turn is 
due to the lack of accurate excited-state potentials. The 
present methodology could have a significant impact in 
this field by providing the necessary information about 
excited-state potentials. 
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Supplementary Online Material Determining 
the Correct Wavefunction out of the Set {'l'i(i)} 

In this supplement we explain how we determine the 
correct wavefunction out of the set of wavefunctions 
i = 1,2,..., 2^, where N is the number of ba- 
sis functions {ipg} needed to span (ref [22] in the 
paper). 

We have defined a set of wavefunctions that can be con- 
structed using the information obtained from the CARS 
signal: 



V'o) 



= Y.\^,)m^{t))^l\^{t)). (14) 



(In writing Eq. (|14p we omitted the proportionality coef- 
ficient ifiEi relative to Eq. (10) in the paper.) Recall that 

^ = Y.g \i'g)(ig \ = Y.g lV's)ag(V'ffl where may take one 
out of two possible values: ±1. A useful property of the 
operator 1 is that its square equals the identity operator 
1: 



11 = '^agag-\^Pg)('lljg\i)g>}{i}. 



a'\ 



gg' 



= ^a2|^,)(V.,|=^|^,)(V,| = l. (15) 



We can derive an equation of motion for ^{t): 

||$(t)) = |l|*(t)) = l||vl/(t))^-*ii/,|*(t)) 

= -ilHjLlYHit)) = -ilHei\^{t)) 

EE -iH,\^{t)), (16) 

where, we have used the fact that 1 is time-independent 
and therefore commutes with Equation shows 
that \E'(i) obeys a time-dependent Schrodinger equation 
with the effective Hamiltonian He = li/el- 

The Hamiltonian He has the conventional form of 
He = Ve + T, where T is the kinetic-energy operator. 
The Hamiltonian Hp therefore takes the form: 



where we emphasize that Ve is time-independent. Let us 
now define the related quantity 



V 



1 



.d_ 



T 



(19) 



where T is the usual kinetic energy operator. Obviously, 
for \&(t) = ^'(t) Eq. ([TO]) is equivalent to the usual time- 
dependent Schrodinger equation for ^'(t) and therefore 
V = Ve is time- independent. We claim that for any other, 
incorrect, wavefunction ^(i), Eq. ([T9l) results in a time- 
dependent potential V. 

In order to show this we substitute T = T + AT in 
Eq. (dl]), where AT = T - T, and obtain: 



Ve = 



1 



= V- 



4-(T + AT) 
[AT]*(i). 



(20) 



The term ^^[AT]^{t) is time- dependent (unless *(i) 
is an eigenf unction of AT, which has no general reason 
to hold. In addition, in the Appendix we show that AT 
is generally different from zero). Therefore, in order to 
preserve the time-independence oi Ve = 1141, V must 
also be time- dependent. 

To summarize: in order to determine the correct wave- 
function out of the set of wavefunctions {'i'i{t)},i = 
1,2,..., 2^, we use the fictitious Schrodinger equation, 
p9)) . to calculate a potential, V, from each wavefunc- 
tion \E'(t) of the set. At different times, t, the wavefunc- 
tions '^{t) will give different potentials V except for the 
one correct wavefunction, ^'(t), that corresponds to the 
correct Schrodinger equation and therefore will give the 
same potential, V = Ve, at all times. Thus, the correct 
wavefunction \['(i) can be selected from the set {^i(i)} 
as the one that provides a time- independent potential via 
Eq. pj)) . Alternatively, as described in the paper, the 
wavefunction ^'(i) that propagates back to the known 
^(0) = ipo, using the corresponding potential calculated 
by Eq. is guaranteed to be the correct reconstructed 
wavefunction, Vl/(i). 



He EE IHet = IVel + ITl = Ve + T , 



(17) 



Note that the operator 1 does not commute with Ve, T 
or He since it does not share a common basis of eigen- 
vectors with the last three operators. Note also that the 
operators Ve, T and He are all time-independent. 
Rearranging Eq. p6)l . we obtain: 



Ve 



1 



'.d ~ 



(18) 



Appendix 

We show that AT ^ 0: 

AT = {f -T) = ll(lTi - T) 

= i(Tl-lT) = l[T,i]. (21) 

The commutator [T, 1] is not identically zero. Therefore, 
1[T, 1] = AT is not identically zero as well. 



